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ABSTRACT 


A  systematic  approach  to  short-period  magnitude  estimation  has  been 
developed  and  applied.  The  approach  uses  newly  developed  statistical 
techniques  in  the  general  linear  model  (GLM)  which  allow  for  the 
problems  of  clipping  and  of  signals  hidden  by  noise.  Measurement 
procedures  are  outlined  and  the  overall  approach  is  first  applied  to 
four  events  in  granite;  PILEDRIVER,  SHOAL,  SAPHIRE,  and  RUBIS.  The 
WWSSN  short-period  network  film  recordings,  with  the  application  of  this 
approach,  form  an  ideal  network  for  shots  over  10  kt  in  hard  rock. 

After  correction  for  the  effects  of  pP ,  estimated  via  syn-hetic  waveform 
calculations,  the  magnitudes  follow  a  theoretical  magni tude : yield  curve 
which  changes  from  a  slope  of  1,0  near  10  kt  to  a  slope  of  0.8  near  100 
kt . 

The  offset  between  the  US  and  Sahara  explosions  is  0.04  !-o  0.12 
units,  with  US  events  biased  low.  It  is  not  clear  from  the  data  if  this 
is  due  to  t*,  or  to  coupling  differences. 

Station  corrections  determined  from  a  suite  of  9  explosions  at 
different  test  sites  around  the  world  show  good  correlation  with 
residuals  estimated  by  North  (1977). 
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Magnitude  (max/GT)  versus  distance  for  PILEJIRIVER.  15 

Downward  and  upward  pointing  arrows  indicate  noise 
and  clipping  limits  respectively.  The  distance 
amplitude  relation  used  is  that  of  Veith  and  Clawson 
(1972)  to  95°  and  of  Sweetser  and  Blandford  (1973) 
beyond,  both  for  zero  kilometers  depth. 

Magnitude  (max/GT)  versus  distance  for  SHOAL.  16 

Downward  and  upward  pointing  arrows  indicate  noise 
and  clipping  limits  respectively.  The  distance 
amplitude  relation  used  is  that  of  Veith  and  Clawson 
(1972)  to  95°  and  of  Sweetser  and  Blandford  (1973) 
beyond,  both  for  zero  kilometers  depth. 

Magnitude  (max/GT)  versus  distance  for  RUBIS.  17 

Downward  and  upward  pointing  arrows  indicate  noise 
and  clipping  limits  respectively.  The  distance 
amplitude  relation  used  is  that  of  Veith  and  Clawson 
(1972)  to  95°  and  of  Sweetser  and  Blandford  (1973) 
beyond,  both  for  zero  kilometers  depth. 

Magnitude  (max/GT)  versus  distance  for  SAPHIRE .  18 

Downward  and  upward  pointing  arrows  indicate  noise 
and  clipping  limits  respectively.  The  distance 
amplitude  relation  used  is  that  of  Veith  and  Clawson 
(1972)  to  95J  and  of  Sweetser  and  Blandford  (1973) 
beyond,  both  for  zero  kilometers  depth. 

Amplitude-distance  relation  for  first  arrivals  as  19 

reported  by  ISC  bulletins  for  January-June  1970. 

Every  event  with  depth  less  than  70  km  and  with 
two  or  more  log^g  (A/T)  values  reported  for 
A  >  110°  were  used.  Occasional  erroneous  first 
picks  which  were  truly  PKP^  were  rejected  on  the 
basis  of  their  time  residuals.  Average  B  factors 
were  computed  by  subtracting  the  log^g  (A/T)  zero- 
to-peak  value  from  the  ISC  m^  value,  averaging  the 
resulting  values  over  2.5°  increments  of  distance, 
and  drawing  a  smooth  curve  by  hand  through  the  re¬ 
sults.  Table  III  gives  the  numerical  B  values  to 
be  used  for  zero  depth,  which  are  obtained  by  add¬ 
ing  0.2.  P,  f  *  an<*  PRIKP  are  indicated  by  small 

d-'.ta.  Initial  arrivals  which  may  be  other  phases 
are  indicated  by  crosses.  The  larger  points  and 
circles  are  the  2.5°  averages.  The  lower  of  the 
double  branch  from  152°— 162®  are  late  picks  of  PKP. 
from  Sweetser  and  Blandford  (1973). 
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WWSSN  synthetic  signals  composed  of  a  direct  23 

and  pP  pulse  from  the  von  Seggern-Blandford 
granite  model  for  the  indicated  yields,  Y, 
delays  and  t*  *  0.4.  For  the  first  column  the 
pP  reflection  coefficient  is  zero,  i.e.  there 
is  no  pP.  For  the  next  two  columns  £he  reflec¬ 
tion  coefficient  is  -0.5  (1  +  exp(-f  ))  where 
f  is  frequency.  The  parameters  T  (period), 

A(log  amplitude),  A/G(log(amplitude/G) )  where 
G  is  the  system  response  at  T,  and  A/GT(log(ampli- 
tude/GT))  are  measured  by  computer  for  the  b  and 
c  phases. 

Experimental  and  theoretical  m.  as  a  function  of  24 

yield  for  the  explosions  PILEDRIVER,  SHOAL,  RUBIS, 
and  SAPHIRE .  The  Worldwide  m^'s  calculated  in 
this  paper  show  an  overall  slope  close  to  1.0. 

When  the  changes  due  to  pP  are  backed  out  (see 
Figure  3  and  Tables  IV,  V  assuming  the  delay 
for  PILEDRIVER  is  0.24  sec)  then  the  slope  be¬ 
comes  even  closer  to  1.0.  The  offset  between 
PILEDRIVER  and  RUBIS  implies  that  the  effects  of 
source  coupling  and  relative  attenuation  between 
the  two  test  sites  result  in  a  bias  of  0.12  m,  with 
the  Sahara  having  the  higher  magnitude  for  a  fixed 
yield.  If  this  offset  is  subtracted  from  the  Sahara 
event  corrected  for  pP  then  the  resulting  amplitude- 
yield  curve  changes  from  a  low-yield  slope  of  1.0 
to  a  high-yield  slope  of  0.8,  in  agreement  with  the 
theoretical  variable  slope,  as  can  be  seen  in  Figure 
6. 


Spectral  ratio  of  the  first  6,4  seconds  of  P  waves  28 
for  PILEDRIVER/SAPHIRE  and  for  PILEDRIVER/RUBIS , 

Note  that  at  NPNT  the  ratio  is  about  30,  whereas  in 
theory  it  should  be  nearly  constant  at  about  0.5, 

This  is  perhaps  due  to  defocussing.  It  is  probably 
not  due  to  differential  Q,  otherwise  there  would  be 
a  greater  slope  to  the  ratio.  For  the  other  ratios, 
the  trend  below  1  Hz  is  toward  the  proper  low-fre¬ 
quency  limit  of  the  yield  ratios  as  indicated  by 
tick  marks  on  the  vertical  axis.  The  change  in  slope 
around  1  Hz  may  reflect  higher  t*  at  NTS  and  poorer 
pP  reflection  at  Hoggar  below  1  Hz  due  to  the  extreme 
relief  0(1:1)  of  the  massif. 
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WWSSN  m.  ' 8  corrected  for  pP  and  source  bias  give  30 

a  variable  slope  ranging  from  1.0  near  SHOAL  to 
0.8  near  SAPHIRE.  For  yields  near  150  kt  the 
correct  slope  is  0.8  after  correction  for  pP. 

This  correct  answer  could  have  been  obtained  in 
several  other  ways,  e.g.  the  simple  slope  between 
RUBIS  and  SAPHIRE  without  correction  for  pP,  and 
(incorrectly)  as  the  simple  slope  from  SHOAL  to 
SAPHIRE  of  the  WWSSN  m. 1 s  averaging  only  detected 
signals . 

Station  corrections  for  the  events  in  Table  VI  33 

(but  using  data  from  the  distance  range  0°<A  <180° 
for  completeness)  as  a  function  of  station  correc¬ 
tions  due  to  North  (1977)  which  were  derived  from 
ISC  earthquake  dates. 
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Film  reading  procedures  used  in  this  study.  12 
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B  factors  for  zero-to-peak  amplitudes  corrected  20 

for  instrument  response  at  period  T  and  divided 

by  T  in  the  distance  range  95°-180°.  To  be  used 

with  maximum  amplitude  in  the  first  5  seconds  of 

the  first  arrival  for  zero  depth  events  (from 

Sweetser  and  Blandford  (1973).  To  be  used  with 

Veith  and  Clawson  (1972)  for  A  <95°. 


25 


m^'s  calculated  using  the  a,  b,  and  max  phase 
for  4  events  for  A  <  95°.  Only  amplitudes  2 
millimeters  on  the  viewing  screen  were  used  as 
signals;  lesser  amplitudes  were  used  as  noise 
estimates,  i.e.,  the  signal  was  judged  to  be 
smaller.  This  is  believed  to  be  the  more  valid 
procedure  due  to  the  difficulty  of  measuring  such 
small  amplitudes,  and  a  tendency  toward  over¬ 
estimation.  If  all  signal  amplitudes  are  used, 
calculations  show  that  the  magnitudes  are  0.02  m^ 
larger  on  the  average. 


Magnitude  summary.  The  column  headed  by  "WW"  is  26 

taken  from  the  last  column  of  Table  IV.  The  next 
column  has  the  pP  effect  removed  but  corrected  so 
that  the  PILEDRIVER  m.  does  not  change.  The  next 

...  D 

column  is  similar  except  that  the  pP  time  is  taken 
to  be  0.17  second  instead  of  0,24  for  PTLEJIRIVER 
[see  Shumway  and  Blandford  (1977)),  and  SHOAL  was 
taken  as  0.15.  (A  special  synthetic  waveform 
calculation  was  performed  for  t=0,15.)  We  note 
a  spread  of  as  much  as  0.08  m^  units  between  these 
last  two  columns.  The  final  column  gives  a  simple 
average  of  signals  neither  below  the  noise  nor 
clipped  on  the  WWSSN, 


n 
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VI  Event  magnitudes  from  entire  WWSSN  network  using  -J2 

max/GT  and  0  <  A  <  95°  with  station  corrections 
common  to  all  events.  The  events  PILEORIVER, 

RUBIS,  SAPHIRE,  and  SHOAL  decrease  in  magnitude 
an  average  of  0.02  units  as  compared  to  Tables 
IV  and  V,  In  this  calculation,  the  absolute 
magnitude  for  these  events  are  more  reliable,  but 
the  relative  magnitudes  less  reliable  since  the 
station  corrections  are  less  finely  tuned  to  the 
NTS  and  Sahara  test  sites.  If  magnitudes  are 
determined  for  0  <  A  <  180°  there  is  typically  an 
increase  by  0,05  m^ ,  except  that  the  Tuomoto  m^  > 
weak  and  strongly  affected  by  PKP,  increase  by 
0,14  mb, 

VII  Station  crustal  corrections  for  WWSSN  stations  36 

based  on  the  crustal  models  given  in  Appendix  II. 

The  corrections  are  the  log  of  the  maximum  peak 
amplitude  of  a  50  kt  explosion  in  granite  with 
a  von  Seggern-B land  ford  source  spectrum,  and  a 
t*  of  0.4  as  seen  through  the  WWSSN  short-period 
response.  The  mean  of  the  logs  has  been  sub¬ 
tracted  out  before  tabulation  so  that  the  net 
effect  on  a  worldwide  network  would  be  zero.  The 
observed  magnitude  should  be  equal  to  the  expected 
magnitude  plus  the  correction. 


INTRODUCTION 


Several  authors  have  approached  the  subject  of  determining  yield 
from  mfe  in  a  systematic  way;  e,g.  Evernden  (  1  967),  Ericsson  (1971a,  b)  , 
Springer  and  Hannon  (1973),  von  Seggern  (1977),  Dahlman  and  Israelson 
(1977),  and  most  recently  Marshall,  Springer,  and  Rodean  (1979). 

Others  have  concentrated  on  the  techniques  of  determining  network 
magnitudes  from  the  individual  station  magnitudes;  among  these  the  work 
of  von  Seggern  (1973),  ~hris tof ferson  et  al.  (1975),  Ringdal  (1976),  and 
von  Seggern  and  Rivers  (1978)  have  been  prominent. 

In  determining  magnitudes,  it  is  most  important  to  have  an  accurate 
distance-amplitude  relationship.  The  earliest  such  relationship  was 
from  Gutenberg  and  Richter  (1956);  in  this  study  we  use  the  relationship 
of  Veith  and  Clawson  (1972),  which  represents  an  average  of  explosion 
data  over  several  test  sites.  When,  on  occasion,  we  choose  to  use  data 
at  distances  larger  than  95  degrees,  we  use  the  relationships  of 
Sweetser  and  Blandford  (1973)  for  PKP  phases. 

The  subject  of  station  effects  at  short-period  stations  has  received 
attention  from  many  authors;  the  most  systematic  collection  of  station 
effects  is  perhaps  that  of  North  (1977),  which  was  obtained  by  analysis 
of  the  readings  reported  to  the  ISC. 

Only  a  few  studies  have  addressed  the  subject  of  how  to  correct 
observed  explosion  magnitudes  for  the  effects  of  pP,  Blandford  et  al. 
(1976)  showed  from  analysis  of  experimental  data  that  the  effect  of  pP 
could  be  large  (  *.3  m.  units),  as  did  von  Seggern  (1977).  Marshall  et 

D 

al,  (1979)  developed  a  systematic  method  of  applying  corrections  for  pP 
and  did  apply  these  corrections  to  a  large  body  of  data. 

In  this  latter  study,  Marshall  et.  al.  also  applied  corrections  for 
differing  levels  of  absorption  under  source  and  receiver.  While  this 
study  was  probably  the  most  careful  and  complete  study  of  magnitude/ 
yield  ever  written,  it  did  suffer  from  a  number  of  weaknesses.  Foremost 
among  these  was  the  small  number  of  observations  for  many  of  the  events. 
This  makes  even  more  serious  the  fact  that  the  study  did  not  use  the 


I 


techniques  of  Ringdai  (1976)  in  order  to  minimize  bias  from  non-detec¬ 
tion  of  small  signals.  Of  course  it  was  impossible  for  Marshall  et  al. 
to  do  this  because  no  mathematical  apparatus  existed  which  made  it 
possible  to  simultaneously  determine  station  corrections  and  take 
account  of  the  biasing  effects  of  noise. 

The  major  advance  of  technique  in  the  present  study  is  to  remove 
this  restriction  which  hampered  Marshall  et  al.,  by  developing  a 
generalization  of  the  general  linear  model  which  can  take  account  of 
noise  and  clipping  levels.  In  addition,  we  actually  model  the  waveforms 
of  the  events  of  interest,  in  order  to  determine  the  pP  correction 
instead  of  simply  using  a  standard  waveform  for  all  yields  as  did 
Marshall  et  al. 


I 

J 
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DATA  MEASUREMENTS 


In  order  to  assure  the  greatest  degree  *  commonality  between 
stations  and  events  in  this  study,  we  have  exclusively  used  the  WWSSN 
short  period  film  data.  While  use  of  this  data  has  been  difficult  in 
the  past  due  to  its  susceptibility  to  either  clipping  or  masking  by 
system  or  earth  noise,  the  advent  of  the  ability  to  allow  for  noise  and 
clipping  within  the  context  of  the  general  linear  model  eliminates  these 
difficulties  and  "transforms"  the  WWSSN  network  into  the  best  possible 
network  for  all  except  small  events.  The  WWSSN  network  becomes  so 
valuable  because  data  with  a  constant  instrument  response  can  be  found 
over  such  a  large  range  of  time  and  with  good  distribution  around  all 
test  sites.  It  is  worth  remarking,  however,  that,  if  possible,  film 
analysis  should  be  confined  to  pre-1978  events  since,  in  that  year,  a 
new  format  of  film  reproduction  began  in  which  there  is  much  less 
resolution. 

The  full  WWSSN  network  was  analyzed  for  each  event,  A  blanket 
request  was  made  for  all  available  data,  and  in  general  terms  for  about 
90%  of  the  stations  film  chips  were  received.  About  20%  of  the  these 
chips  were  unusable  because  of  poor  film  development,  timing  failure, 
missing  7,  component,  etc. 

Table  I  shows  the  film  reading  procedures  that  were  given  to  the 
analysts  to  ensure  that  proper  noise  readings  were  made  if  the  signal 
was  not  visible. 


TABLE  I 


Film  Reading  Procedures  Used  in  this  Study 

Read  the  a  (zero-crossing  to  first  peak),  b  (first  peak  to  first 
trough),  and  max  peak-to-trough  or  t rough-to-peak  in  first  5  seconds 
in  millimeters.  Record  readings  in  millimeters  in  the  data  file. 
Also,  record  gain  as  written  on  film  and  as  seen  with  particular 
viewer . 

Measure  period  for  max  as  peak-to-peak  or  zero-crossing  to  zero¬ 
crossing  or  trough-to-trough  as  in  your  judgement  best  reflects 
the  period  of  the  maximum  energy. 

If  a  weak  signal  must  be  measured  then  try  to  find  a  strong 
signal  from  the  same  test  site  and  by  correlation  try  to  establish 
if  a  particular  cycle  is  a  or  b.  In  general,  it  is  good  practice 
to  analyze  events  in  pairs,  with  any  weak  events  paired  with  a 
strong  event. 

If  it  is  not  possible  to  determine  the  exact  location  of  the  a  or 
b  cycles  then  the  first  clear  down-swing  is  to  be  used  to  place  a 
"noise"  limit  on  the  b  cycle;  that  is,  the  true  b  amplitude  is  less 
than  or  equal  to  that  down-swing.  Similarly  for  the  a  phase  if 
there  is  a  clear  up-swing  which  cannot  confidently  be  said  to  be 
the  a  phase. 

More  commonly,  however,  there  will  be  no  clear  up-swing;  in  this 
case  a  noise  measurement  must  be  made,  and  to  do  so  search  the 
preceding  20  seconds  for  the  largest  peak-to-trough  excursion  in  the 
1-2  Hz  frequency  range.  (For  the  a  phase  find  the  largest  zero-to- 
peak.)  Often  such  amplitudes  are  less  than  1  mm  on  WWSSN  film.  The 
period  is  to  be  recorded  as  1  Hz. 

Of  course,  if  no  arrival  at  all  can  be  discerned  then  again  a  noise 

measurement  must  be  made. 

If  the  data  for  any  phase  are  "wiped  out"  then  an  estimate  of  the 
clipping  level  must  be  made.  To  do  this  measure  the  amplitude  of 
the  largest  "turning  point"  visible  and  multiply  by  2.  If  you 
can  be  confident  that  the  largest  turning  point  is  off  the  film 
then  you  could  use  the  maximum  distance  to  the  edge  of  the  film. 


DATA  ANALYSIS 


In  Table  II  we  see  the  events  studied  in  detail  in  this  report,  and 
in  Figures  la-d  we  see  the  magnitude  versus  distance  plots  for  those  ' 

events.  The  upward  directed  arrows  indicate  that  the  signal  clipped  and 
that  the  true  magnitude  is  greater  than  that  value;  a  downward  directed 
arrow  indicates  a  noise  reading  and  that  the  magnitude  is  less  than  that 
value.  We  see  that  the  data  are  consistent  with  the  hypothesis  that  the 
amplitude  distance  curves  are  correct  and  that  the  magnitude  is 
statistically  constant  as  a  function  of  distance. 

Figure  2  is  a  plot  of  the  amplitude-distance  relationship  for 
shallow  earthquakes  at  typical  depths  of  40  kilometers  and  for  distances 
greater  than  90  degrees;  in  Table  III  the  actual  B  factors  used  for  0  km 
depth  are  tabulated.  These  are  to  be  used  together  with  the  Veith  and 
Clawson  curves  for  0  depth.  (Table  III  is  for  zero-to-peak  measure¬ 
ments;  it  must  be  remembered  that  Veith  and  Clawson's  tables  are  for 
peak-to-peak. ) 

The  general  linear  model  used  on  the  measured  noise  and  signal 
amplitude  data  is  discussed  in  detail  in  Appendix  I.  Here  we  mention 
only  that  the  method  is  an  iterative  one  in  which  initial  guesses  as  to 
the  event  magnitudes,  station  corrections,  and  magnitude  standard  devia¬ 
tions  are  updated  in  such  a  way  that  the  likelihood  of  the  new  estimates 
is  monotonies  1 ly  increasing.  Also,  it  is  worth  mentioning  that  at  each 
stage  in  the  iteration,  the  amplitudes  at  those  stations  that  do  not 
observe  the  signal  are  estimated  as  being  equal  to  the  mean  value  so  far 
computed,  corrected  by  the  station  correction  as  so  far  computed,  and 
adjusted  by  the  most  probable  level  of  the  true  observation  below  noise 
or  above  the  clipping  level.  These  amplitudes  are  averaged  at  each 
stage  in  the  iteration  in  order  to  give  the  magnitude  estimate  for  the 
next  stage  in  the  iteration.  The  principal  reference  for  this  new 
procedure  is  that  of  Dempster  (1977)  and  the  extension  to  clipped  and 
noise-hidden  data  was  done  by  Aitkin  (1981)  and  by  Schmee  and  Hahn 
(1979). 

?! 
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TABLE  II 


Events  Analyzed  in  this  Study 


EVENT 

DATE 

TIME 

LOCATION 

DEPTH  (m) 

YIELD 

PILEDRIVER 

66/06/02 

15:30:00 

37 . 07  °N 

116. 07°W 

462 

62 

SHOAL 

63/10/26 

17:00:00 

39 . 20°N 

118. 30°E 

367 

12 

RUBIS 

63/10/20 

13:00:00 

24:02:07. 8°N 
05:02: 19. 0°E 

— 

52 

SAPHIRE 

65/02/27 

11/30:00 

24:03:31. 4°N 
05:01:52. 3°E 

645-785 

120 

Yields  from  Marshall  et  al  (1979),  depth  of  SAPHIRE  from  I^EA 
publication,  November  1970.  Locations  of  French  shots  from  Bolt 
(1976). 


MAGNITUDE 


Figure  la.  Magnitude  (max/GT)  versus  distance  for  PILEDRIVER. 

Downward  and  upward  pointing  arrows  indicate  noise 
and  clipping  limits  respectively.  The  distance 
amplitude  relation  used  is  that  of  Veith  and  Clawson 
(1972)  to  95°  and  of  Sweetser  and  Blandford  (1973) 
beyond,  both  for  zero  kilometers  depth. 
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Figure  lb.  Magnitude  (max/GT)  versus  distance  for  SHOAL.  Downward  and 
upward  pointing  arrows  indicate  noise  and  clipping  limits 
respectively.  The  distance  amplitude  relation  used  is  that  of 
Veith  and  Clawson  (1972)  to  95®  and  of  Sweetser  and  Blandford 
(1973)  beyond,  both  for  zero  kilometers  depth. 
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MAGNITUDE 
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Figure  lc.  Magnitude  (max/GT)  versus  distance  for  R1JBIS. 

Downward  and  upward  pointing  arrows  indicate  noise 
and  clipping  limits  respectively.  The  distance 
amplitude  relation  used  is  that  of  Veith  and  Clawson 
(1972)  to  95"  and  of  Sweetser  and  Blandford  (1973) 
beyond,  both  for  zero  kilometers  depth. 
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Figure  Id,  Magnitude  (max/GT)  versus  distance  for  SAPHIRE. 

Downward  and  upward  pointing  arrows  indicate  noise 
and  clipping  limits  respectively.  The  distance 
amplitude  relation  used  is  that  of  Veith  and  Clawson 
(1972)  to  95°  and  of  Sweetser  and  Blandford  (1973) 
beyond,  both  for  zero  kilometers  depth. 
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TABLE  III 


B  factors 

for 

zero-to-peak 

amplitudes 

corrected 

for  instrument  response  at 

period  T  and  divided 

by  T  in  the  distance  range 

95°-180° . 

To  be 

used 

with  maximum 

amplitude  in 

the  first  5 

seconds  of 

the  first 

arrival  for  zero 

depth  events  (from 

Sweet ser 

and 

Blandford  (1973).  To  be 

used 

with 

Veith  and 

Clawson  (1972)  f 

or  A  <  95°  . 

Distance 

Distance 

(deg) 

B 

(deg) 

B 

95 

4.430 

137 

4.600 

96 

4.500 

138 

4.590 

97 

4.560 

139 

4.575 

98 

4.625 

140 

4.565 

99 

4.680 

141 

4.550 

100 

4.735 

142 

4.450 

101 

4.790 

143 

4.215 

102 

4.840 

144 

3.850 

103 

4.905 

145 

3.710 

104 

4.980 

146 

3.490 

105 

5.060 

147 

3.675 

106 

5.180 

148 

3.680 

107 

5.325 

149 

3.700 

108 

5.425 

150 

3.725 

109 

5.400 

151 

3.750 

110 

5.225 

152 

4.180 

111 

5.020 

153 

4.210 

112 

4.950 

154 

4.235 

113 

4.915 

155 

4.260 

114 

4.875 

156 

4.275 

115 

4.850 

157 

4.300 

116 

4.830 

158 

4.340 

117 

4.800 

159 

4.370 

118 

4.770 

160 

4.375 

119 

4.730 

161 

4.385 

120 

4.690 

162 

4.385 

121 

4.660 

163 

4.385 

122 

4.625 

164 

4.375 

123 

4.585 

165 

4.375 

124 

4.550 

166 

4.375 

125 

4.515 

167 

4.375 

126 

4.480 

168 

4.370 

127 

4.450 

169 

4.360 

128 

4.425 

170 

4.360 

129 

4.415 

171 

4.360 

130 

4.425 

172 

4.360 

131 

4.450 

173 

4.360 

132 

4.490 

174 

4.360 

133 

4.550 

175 

4.360 

134 

4.590 

176 

4.360 

135 

4.610 

177 

4.360 

136 

4.610 

178 

4.360 
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Application  of  this  analysis  procedure  to  the  measurements  for  the 
four  events  in  Table  II  yields  the  magnitudes  in  Table  IV.  Here  we  see 
that  we  have  imposed  the  limitation  that  the  distance  be  less  than  95 
degrees  in  order  to  conform  to  the  standard  definition  of  magnitude. 
We  also  see  that  magnitudes  have  been  calculated  for  the  a,  b,  and  max 
phase,  using  only  the  raw  amplitudes,  the  amplitudes  corrected  for 
instrument  response  at  the  measured  period,  and  corrected  for  instrument 
response  and  divided  by  period.  For  the  max  phase,  this  last 
measurement  is  in  accord  with  the  standard  definition  of  and  it  is 
this  magnitude  only  which  we  shall  discuss  in  the  rest  of  this  report, 
although  for  all  events,  the  measurements  have  been  made  to  enable 
calculation  of  the  other  magnitudes.  The  relationship  between  the 
magnitudes  remains  a  subject  for  further  research. 

We  may  remark,  however,  on  a  few  points  of  interest  which  can  be 
seen  on  inspection  of  Table  IV. 

o  Moving  from  the  a  to  b  phase,  there  is  an  increase  of 
about  0.6  m^  units;  from  the  b  to  max  phase,  there  is  a 
difference  of  about  0.25  m^  units,  for  a  total  difference 
from  a  to  max  of  about  0.85  units. 

o  For  all  events  except  SHOAL,  the  magnitude  increases  slightly 
when  corrected  for  instrument  response,  and  then  decreases 
when  divided  by  period.  For  SHOAL  it  is  the  reverse,  which 
may  have  to  do  with  the  dominant  period  of  this  smaller 
event.  For  SHOAL  the  magnitude,  when  corrected  for  in¬ 
strument  response  and  divided  by  period,  is  nearly  identical 
to  the  case  of  magnitudes  computed  from  raw  amplitudes. 
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CORRECTION  FOR  t*  AND  FOR  pP  AND  OTHER  SOURCE  EFFECTS 

Before  the  m^  values  reported  in  Table  IV  can  be  accurately  related 
to  yield,  it  is  necessary  to  make  corrections  for  the  effects  of  pP  and 
of  absorption,  much  as  was  done  by  Marshall  et  al,  (1979).  To  perform 
calculations  of  synthetic  waveforms  to  evaluate  the  effects  of  pP,  it  is 
necessary  to  have  the  delay  times.  The  delay  times  for  RUBIS  and 
SAPHIRE  (  t=0.21,  0,27  sec)  can  be  obtained  from  a  1970  IAEA 
publication,  which  also  gives  the  velocity  of  the  French  granite  as  ^.25 
km/ sec  (P.  Marshall,  personal  communication),  Marshall  et  al,  (1979) 
also  give  the  delay  times  for  PILEDRIVER  and  SHOAL  as  (0.24,  0.21 
seconds);  however,  Shumway  and  Blandford  (1977)  found  the  delay  for 
PILEDRIVER  to  be  0,17  seconds;  and  cube  root  depth  scaling  would  give 
the  delay  for  SHOAL  to  be  0,15  seconds.  Both  sets  of  delays  will  be 
used  in  this  study,  and  we  shall  see  that  there  are  significant  effects 
resulting  from  such  seemingly  minor  differences  in  pP  delays. 

Typical  results  of  such  calculations  may  be  seen  in  Figure  5;  the  pP 
corrections  may  be  obtained  by  subtracting  the  magnitude  in  the  first 
column  from  the  appropriate  magnitude  seen  in  the  second  or  third 
column.  The  results  of  such  calculations  are  seen  in  Table  V.  We  have 
chosen  to  correct  all  of  the  events  for  pP  and  then  to  add  a  constant  so 
that  the  PILEDRIVER  m^  remains  unchanged. 

For  example,  the  PILEDRIVER  magnitude  for  the  c  phase  corrected  for 
instrument  response  and  divided  by  T  with  the  pP  phase  included  is,  from 
Figure  3,  3.65  (to  within  an  arbitrary  constant).  Without  pP  the 
magnitude  is  3.42,  The  difference  is  a  decrease  of  0.23  magnitude 
units.  Now  SAPHIRE  decreases  from  3,84  to  3,64,  a  decrease  of  only  0.19 
units  so  that  the  0.23  correction  applied  to  SAPHIRE  results  in  an 
increase  from  5.79  in  the  fourth  column  of  figures  to  5.83  in  the  fifth. 
Notice  also  the  significant  differences  between  columns  5  and  6,  for 
example  changes  of  .08  m.  for  RUBIS  and  SAPHIRE.  This  shows  that 

D 

apparently  minor  differences  in  pP  delay  times  can  have  significant 

effects  on  m,  . 

D 

In  Figure  4  we  see  plotted  some  of  the  m,  values  from  Table  V.  We 
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Figure  4.  Experimental  and  theoretical  m.  as  a  function  of  yield  for  the 

explosions  PILEDRIVER,  SHOAL,  RUBIS,  and  SAPHIRE.  The  Worldwide 
m^'s  calculated  in  this  paper  show  an  overall  slope  close  to 
1.0.  When  the  changes  due  to  pP  are  backed  out  (see  Figure  3 
and  Tables  IV,  V  assuming  the  delay  for  PILEDRIVER  is  0,24  sec) 
then  the  slope  becomes  even  closer  to  1,0.  The  offset  between 
PILEDRIVER  and  RUBIS  implies  that  the  effects  of  source  coupling 
and  relative  attenuation  between  the  two  test  sites  result  in  a 
bias  of  0.12  m^  with  the  Sahara  having  the  higher  magnitude  for 
a  fixed  yield.  If  this  offset  is  subtracted  from  the  Sahara 
event  corrected  for  pP  then  the  resulting  amplitude-yield  curve 
changes  from  a  low-yield  slope  of  1.0  to  a  high-yield  slope  of 
0.8,  in  agreement  with  the  theoretical  variable  slope,  as  can  be 
seen  in  Figure  6, 
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see  that  the  simple  values  (m^  values  which  are  the  simple  average  of 
the  detecting  station  magnitudes)  have  a  slope  of  close  to  0.8  as  a 
function  of  yield.  This  effect  is  mostly  due  to  the  overestimate  of  the 
SHOAL  magnitude.  When  the  estimates  which  take  account  of  noise  and 
clipping  levels  are  included,  however,  the  slope  becomes  close  to  1.0; 
and  when  the  further  corrections  are  made  for  pP ,  assuming  the  pP  delay 
for  PILEDRIVER  is  0.24  seconds,  then  the  slight  increase  in  m^  for 
SAPHIRE  leads  to  an  even  closer  fit  to  a  1,0  slope. 

The  dashed  lines  touching  the  RUBIS  and  PILEDRIVER  events  indicate 
an  offset  of  0.12  m^  units  between  the  NTS  and  Sahara  test  sites.  (Use 
of  a  pP  delay  for  PILEDRIVER  of  0.17  seconds  would  yield  an  offset  of 
only  0.04  m^ , )  We  may  ask  if  this  offset  could  be  due  to  absorption. 

In  Figure  5  we  see  the  ratios  of  PILEDRIVER  spectra  divided  by 
spectra  due  to  either  SAPHIRE  or  RUBIS,  Considering  the  spectra  as  a 
whole,  one  could  not  reject  the  hypothesis  that  the  trend  with  frequency 
was  due  to  a  relative  t*  between  the  test  sites  of  t*  =  0,06;  a  value 
which  would  explain  an  offset  of  0,08  m.  units.  However,  it  is  notable 

D 

that  between  about  1  and  2.5  Hz  there  is  quite  a  steep  slope 

corresponding  to  a  relative  t*  of  about  0.4.  If  this  truly  represented 

the  difference  in  t*  between  NTS  and  the  Sahara,  then  one  would  expect 

the  NTS  m.  to  be  as  much  as  0,6  m.  lower  than  the  Sahara  m.  instead  of 
D  D  u 

the  observed  0,12  m^,  Thus  one  is  immediately  tempted  to  assert  that 
the  apparent  large  slope  between  1  and  2.5  Hz  is  simply  an  artifact. 

On  the  other  hand,  DsMman  and  Israelson  (  1977)  write  that  the 
cavity  size  from  the  French  tests  was  substantially  smaller  than  that 
expected  on  the  basis  of  NTS  tests.  If  so,  then  the  French  tests  may  be 
small  because  of  some  unknown  coupling  factor;  and  the  NTS  explosions 
coincidentally  small  because  of  high  absorption.  The  data  basically  are 
insufficient  to  resolve  such  complicated  questions. 

As  a  marginal  remark,  one  may  speculate  that  if  there  were  indeed 
substantially  larger  attenuation  under  NTS  than  under  the  Sahara,  then 
the  turnover  of  the  spectral  ratio  near  1  Hz  in  Figure  5  may  be  due  to 
the  poor  reflection  of  pP  in  the  Sahara  due  to  the  fact  that  the 
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AMPLITUDE  RATIO 


FREQUENCY  Hi 

0.0  1.0  2.0  3.0  4.0  5.0 
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Figure  5.  Spectral  ratio  of  the  first  6.4  seconds  of  P  waves  for 

PILEDRIVER/SAPHIRE  and  for  PILEDRIVER/RUBIS .  Note  that  at  NPNT 
the  ratio  is  about  30,  whereas  in  theory  it  should  be  nearly 
constant  at  about  0,5.  This  is  perhaps  due  to  defocussing.  It 
is  probably  not  due  to  differential  Q,  otherwise  there  would  be 
a  greater  slope  to  the  ratio.  For  the  other  ratios,  the  trend 
below  1  Hz  is  toward  the  proper  low-frequency  limit  of  the  yield 
ratios  as  indicated  by  tick  marks  on  the  vertical  axis.  The 
change  in  slope  around  1  Hz  may  reflect  higher  t*  at  NTS  and 
poorer  pP  reflection  at  Hoggar  below  1  Hz  due  to  the  extreme 
relief  0(1:1)  of  the  massif. 
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explosions  were  detonated  in  the  Hoggar  massif,  which  in  profile  is 
approximately  1  km  high  and  2  kilometers  wide.  If  the  reflection  is 
poor,  then  the  P-wave  amplitude  at  low  frequencies  from  the  Hoggar  would 
not  be  efficiently  cancelled  by  pP  and  the  NTS/Sahara  ratio  would 
decline  toward  low  frequencies,  as  observed.  If  this  were  the  case  than 
one  would  deduce  that  poor  coupling  in  the  Sahara  is  being  offset  by 
high  absoption  under  NTS. 

Returning  to  Figure  4,  we  let  the  source  of  the  0.12  m^  offset 
remain  undecided  and  simply  shift  the  two  Sahara  shots  down  by  0.12  m^ . 
The  result  is  seen  in  Figure  6  and  we  see  that  the  observed  m^tyield 
curve  is  similar  in  shape  to  the  theoretical  curve.  The  only  major 
discrepancy  is  that  the  SHOAL  event  m^  is  too  small.  This  may,  perhaps, 
be  within  the  margin  of  error  for  this  small  event  at  a  site  rather 


Figure  6.  WWSSN  m  's  corrected  for  pP  and  source  bias  give 
a  variable  slope  ranging  from  1.0  near  SHOAL  to 
0.8  near  SAPHIRE.  For  yields  near  150  kt  the 
correct  slope  is  0.8  after  correction  for  pP , 

This  correct  answer  could  have  been  obtained  in 
several  other  ways,  e.g.  the  simple  slope  between 
RUBIS  and  SAPHIRE  without  correction  for  pP,  and 
(incorrectly)  as  the  simple  slope  from  SHOAL  to 
SAPHIRE  of  the  WWSSN  *  s  averaging  only  detected 
signals. 
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STATION  EFFECTS 


The  station  effects  which  were  calculated  in  the  course  of  producing 
the  magnitudes  in  Table  V  and  Figures  4  and  6  are  the  appropriate  ones 
for  a  comparison  of  magnitudes  between  the  NTS  and  Sahara  test  sites. 

If  one  wanted  to  add  one  more  test  site,  e,g,,  Semipalat insk ,  then  there 
would  be  a  different  set  of  station  corrections.  However,  we  may  be 
interested  in  an  estimate  of  station  corrections  for  which  the 
peculiarities  of  the  source-station  paths  have  been  averaged  out, 
resulting  in  station  corrections  which  fairly  closely  represent  the 
near-receiver  structure  and  absorption.  For  this  purpose  one  would  want 
to  have  as  well-distributed  a  set  of  events  as  possible,  without  overly 
weighting  any  particular  test  site. 

For  this  purpose,  we  have  measured  r entire  WWSSN  network  for  the 
events  listed  in  Table  VI,  The  first  9  of  these  events  were  then  used 
to  determine  station  corrections,  and  these  station  corrections  were 
used  to  determine  the  magnitudes  of  the  last  five  events  in  Table  VI. 
Actually  the  Table  VI  magnitudes  and  corresponding  station  corrections 
were  determined  using  only  readings  for  distances  less  than  95  degrees. 
Thus  these  magnitudes  are  perhaps  the  most  authoritative  available  in  an 
absolute  sense.  We  note  that  there  are  small  changes,  averaging  a 
decrease  of  0.02  m^  units  as  compared  to  the  magnitudes  in  Table  V,  It 
is  important  to  remember,  however,  that  the  magnitudes  in  Table  V  are 
the  ones  to  use  for  comparison  of  the  NTS  and  Sahara  test  sites. 

To  determine  station  corrections,  however,  it  is  useful  to  allow  PKP 
readings;  otherwise  many  stations  do  not  have  sufficient  readings  to 
define  the  station  effects.  Carrying  this  through  for  the  first  9 
events  in  Table  VI,  we  plot  the  results  versus  the  station  corrections 
determined  by  North  (1977),  as  seen  in  Figure  7.  By  no  means  all  of  the 
WWSSN  stations  can  be  found  in  North's  list,  and  North  determined 
station  corrections  for  many  non-WWSSN  stations.  In  Figure  7  we  see 
that,  in  general,  there  is  a  definite  correlation;  however,  there  are 
several  stations  for  which  the  North  corrections  are  too  large.  We 
interpret  this  to  show  that  the  poorer  stations  tend  to  detect  only  the 
anomalously  large  amplitudes.  As  a  result,  North's  analysis,  which 
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TABLE  VI 


Event  magnitudes  from  entire  WWSSN  network  using  max/GT  and  0  <  A  < 95°  with 
station  corrections  common  to  all  events.  The  events  PTLEDRIVER,  RUBIS, 
SAPHIRE,  and  SHOAL  decrease  in  magnitude  an  average  of  0.02  units  as 
compared  to  Tables  IV  and  V,  In  this  calculation,  the  absolute  magnitude 
for  these  events  are  more  reliable,  but  the  relative  magnitudes  less 
reliable  since  the  station  corrections  are  less  finely  tuned  to  '  e  NTS  and 
Sahara  test  sites.  If  magnitudes  are  determined  for  0 <  A  <  180°  there  is 

typically  an  increase  by  0.05  m^ ,  except  that  the  Tuomoto  m^ ,  weak  and 
strongly  affected  by  PKP,  increase  by  0.14  m^ , 


Event 

Oa  te 

Latitude 

Longi tude 

mb 

Piledriver 

66/06/02 

37 . 23aH 

1 16 , 06°  W 

5.49 

Rubis 

63/10/20 

24 . 00°  N 

5 . 00°E 

5.45 

Azgir 

76/07/29 

47 . 782°N 

48. 1  20°  E 

5.84 

Novaya  Zemlya 

76/10/20 

7  3 , 00°  N 

5  5 . 00°  E 

4.76 

Kazakh 

73/12/14 

50 . 04°  N 

79 . 01  “E 

5.76 

Tuomoto 

77/02/19 

22 . 1 0°N 

13«. 76°W 

4.66 

Ch  ina 

76/10/17 

41 ,00°N 

89 , 00°  E 

4.56 

Sa  lmon 

64/10/22 

31 .  14°N 

89. 57°W 

4.40 

Longshot 

65/10/29 

51 ,44°N 

179. 18°E 

5 .  SO 

Shoa  1 

63/10/26 

39. 20°N 

1 1 8 . 38° W 

4.75 

Saphire 

65/02/27 

24 , 06°  N 

5 . 03°E 

5.7? 

Kazakh 

69/09/11 

49. 70“N 

78. 1 1°E 

4.57 

Kazakh 

78/09/15 

49 .91 °N 

78 . 94°E 

5.76 

Azgir 

77/09/30 

4  7 . 800°  N 

48. 145°E 

4.81 

Station  corrections  determined 

with  first  9 

events  in  list. 
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(max/GT)  STATION  CORRECTIONS 


could  make  no  allowance  for  noise,  would  tend  to  overestimate  the 
station  correction.  Some  support  for  this  may  be  garnered  from  the  fact 
that  stations  which  are  generally  regarded  as  sensitive  and  well-run, 
such  as  ALQ,  GOL,  COL,  KBL,  BKS,  and  KON  do  tend  to  lie  close  to  the 
line  through  the  origin  with  slope  1.0. 

Improvement  of  these  station  corrections  could  be  achieved  by  such 
approaches  as  adding  data  from  Gasbuggy,  Rulison,  Faultless,  a  bigger 
event  from  Tuamoto,  Soviet  PNE's,  and  shots  from  Novaya  Zemlya,  It 
would  be  useful  also  to  add  more  events  from  the  established  test  sites 
in  order  to  fill  in  stations  which  did  not  observe  the  selected  event 
and  to  average  out  fluctuations  which  come  even  from  within  a  test  site. 
However,  an  extension  of  the  present  program  is  needed  before  this  last 
calculation  can  be  performed;  for  each  event  from  a  particular  test  site 
weights  roughly  inversely  proportional  to  the  number  of  events  from  that 
test  site  must  be  assigned  in  order  to  prevent  overweighting  of  station 
corrections  by  a  particular  test  site.  This  apparently  simple 
modification  calls  for  rather  substantial  changes  in  the  computation 
techniques , 

In  a  preliminary  attempt  to  see  at  what  distance  between  events  the 
station  corrections  became  independent  of  each  other,  station 
corrections  were  computed  for  pairs  of  events  at  increasing  distances 
from  each  other.  As  the  distance  between  events  increased,  it  was 
expected  that  the  residual  variance  would  increase.  We  found  that  this 
was  the  case,  and  also  that  it  was  necessary,  in  order  to  obtain  an 
unbiased  estimate  of  the  standard  deviation,  to  divide  the  total 
variance  not  by  the  number  of  stations  x  events,  but  by  the  number  of 
observations,  including  noise  and  clipping  observations. 

When  this  was  done,  we  found  that  for  the  two  event  pairs  with 
separations  of  approximately  3  kilometers  (RUBIS-SAPHIR  and  AZGIR 
07/76-09/77),  the  standard  deviation  was  about  0,12  m^;  but  that  for 
event  pair  separations  from  15  kilometers  out  to  100  degrees,  the 
computed  residual  standard  deviations  were  not  statistically  different 
from  0.3  m^ . 


As  a  final  task  in  this  study  we  surveyed  the  geological  literature 
to  find  crustal  models  for  each  of  the  WWSSN  stations.  The  details  of 
these  crustal  models  may  be  found  in  Appendix  II.  Using  these  crustal 
models,  a  synthetic  waveform  was  computed  using  a  causal  wavelet  from  a 
50  kt  explosion  through  a  t*  of  0.4  and  a  WWSSN  SP  response.  The  peak 
to  trough  maximum  amplitude  of  the  waveforms  were  measured,  the 
logarithms  were  taken  and  the  mean  removed.  The  results  are  displayed 
in  Table  VII.  Application  of  these  corrections  should  reduce  the 
variance  of  the  magnitude  residuals,  but  will  not  change  the  mean  value 
if  observations  are  made  at  all  stations.  Future  research  should  test 
whether  or  not  these  corrections  do  in  fact  reduce  the  observed 
variance . 
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TABLE  VII 


Station  crustal  corrections  for  WWSSN  stations  based  on  the  crustal  models 
given  in  Appendix  II,  The  corrections  are  the  log  of  the  maximum  peak 
amplitude  of  a  50  kt  explosion  in  granite  with  a  von  Seggern-Bland ford 
source  spectrum,  and  a  t*  of  0,4  as  seen  through  the  WWSSN  short-period 
response.  The  mean  of  the  logs  has  been  subtracted  out  before  tabulation 
so  that  the  net  effect  on  a  worldwide  network  would  be  zero.  The  obstrved 
magnitude  should  be  equal  to  the  expected  magnitude  plus  the  correction. 


AAE 

-0.061 

DAV 

-0.048 

LPB 

0.071 

RCD 

0.160 

A  AM 

0.069 

DUG 

-0.059 

LPS 

-0.059 

RIV 

0.019 

ADE 

-0.059 

ESK 

0.001 

LUB 

0.079 

RKO 

-0.059 

AFI 

-0.059 

FLO 

-0.040 

MAL 

-0.020 

SBA 

-0.059 

AKU 

-0.059 

GDH 

-0.059 

MAN 

-0.059 

SCP 

-0.059 

ALQ 

-0.059 

GEO 

-0.059 

MAT 

-0.048 

SDB 

-0.059 

ANP 

-0.047 

GIE 

-0.048 

MDS 

-0.059 

SEO 

-0.059 

ANT 

-0.020 

GOL 

-0.059 

MNN 

-0.022 

SHA 

0.247 

AQU 

-0.036 

GSC 

-0.0  59 

MSH 

-0.059 

SHI 

0.060 

ARE 

-0.040 

GUA 

-0.037 

MUN 

-0.059 

SHK 

-0.059 

ATL 

-0.059 

HKC 

-0.0  59 

NAI 

-0.059 

SHL 

0.001 

ATU 

0.001 

HLW 

0.056 

NAT 

0.079 

SJG 

-0.048 

BAG 

-0.023 

HNM 

-0.036 

NDI 

0.059 

SNA 

0.252 

BEC 

0.142 

HNR 

0.050 

NNA 

-0.059 

SNG 

-0.059 

BHP 

0.140 

HOW 

0.124 

NOR 

0.050 

SOM 

0.079 

BKS 

-0.008 

1ST 

-0.059 

NUR 

-0.059 

SPA 

0.100 

BLA 

-0.148 

JCT 

0.056 

0GD 

-0.058 

STU 

0.001 

BOG 

0.041 

JER 

0.079 

0XF 

0.079 

TAB 

-0.048 

BOZ 

0.054 

KEV 

-0.059 

PDA 

0.026 

TAU 

0.116 

BUL 

-0 .007 

KIP 

-0.059 

PEL 

0.050 

TOL 

0.208 

CAR 

0.050 

KOD 

0.124 

PMG 

0.249 

TRI 

0.079 

CHG 

-0.059 

KON 

-0.059 

POO 

-0.059 

TRN 

-0.059 

CMC 

-0.059 

KTG 

-0.059 

PRE 

-0.059 

TUC 

0.046 

COL 

-0.059 

LAH 

0.259 

PTO 

-0.059 

UME 

-0.059 

COP 

-0.016 

LEM 

-0.048 

QUE 

0.050 

VAL 

-0.059 

COR 

0.075 

LON 

-0.059 

QUI 

-0.065 

WEL 

0.0  56 

CTA 

-0.059 

LOR 

0.050 

RAB 

-0.059 

WES 

-0.0  59 

DAL 

-0.018 

LPA 

0.124 

RAR 

-0.048 

WIN 

-0.059 
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SUMMARY  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 


By  use  of  newly  developed  mathematical  analysis  techniques,  the  full 
WWSSN  short-period  network  has  been  used  to  obtain  unbiased  magnitude 
estimates  for  four  shots  in  granite;  PILEDRIVER,  SHOAL,  SAPHIRF. ,  and 
RUBIS.  After  correction  for  pP,  there  appears  to  be  an  offset  between 
the  US  and  Sahara  shots  of  0.04  to  0.12  magnitude  units.  If  this  offset 
is  eliminated,  then  the  magni tude / y ie  Id  curve  agrees  with  a  theoretical 
curve  in  that  the  slope  changes  from  1.0  from  10  to  50  kt,  to  0.8  from 
50  to  100  kt.  The  data  are  inadequate  to  determine  if  this  offset  is  a 
result  of  mantle  absorption,  poor  pP  reflection  at  low  frequencies,  or 
differences  in  coupling  between  US  and  Sahara  granite. 

Magnitudes  have  also  been  estimated  for  10  additional  events 
including  events  from  Azgir,  Novaya  Zemlya,  East  Kazakh,  Tuomoto,  China, 
Mississippi  (Salmon)  and  Amchitka  (Longshot),  The  station  corrections 
have  been  determined  s imul taneous ly  for  all  these  source  mechanisms  so 
that  the  magnitudes  of  all  of  these  source  regions  are  connected  to  each 
other  in  a  consistent  manner. 

The  station  corrections  determined  from  this  suite  of  4  source 
regions  have  been  plotted  versus  the  corrections  of  North  (1977)  and 
show  substantial  correlation.  The  correlation  is  best  for  those 
stations  thought  to  be  the  best  run.  This  is  in  agreement  with  the  idea 
that  at  the  poorer  stations  only  the  abnormally  large  arrivals  would  be 
reported  to  the  ISC,  the  source  of  the  data  used  to  determine  North’s 
cor rec t ions , 

Further  research  should  proceed  along  the  lines  of: 

o  Determine  pP  delay  times  directly  from  the  data  instead 
of  from  a  priori  information. 

o  Generalize  the  mathematical  procedures  so  that  several 
events  can  be  used  from  a  test  site  without  biasing  the 
station  corrections. 
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o  Apply  a  priori  station  corrections  and  see  if  the  estimated 
variance  decreases. 

o  Better  delineate  the  distance  at  which  the  residual 

variance  after  station  corrections  increases  from  0.12 

to  0.3  m.  units, 
b 

o  Further  check  the  idea  that  the  station  corrections  de¬ 
termined  in  this  study  are  superior  to  those  from  earlier 
studies  by  using  an  objective  measure  of  station  quality 
such  as  reporting  threshold  and  plotting  this  measure  versus 
the  difference  between  the  corrections  determined  in  this 
study  and  those  determined  in  other  studies. 

o  Measure  additional  events  with  known  yields  such  as  Gasbuggy, 
Faultless,  and  Rulison;  and  improve  the  estimation  of  station 
effects  by  measuring  shots  in  other  test  sites. 

o  Extend  the  analysis  to  Mg  magnitude. 
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Joint  Maximum  Likelihood  Estimation  of  Seismic 
Magnitude  and  Distance-Amplitude  Dependence 
In  the  Presence  of  Clipped  and  Missed  Signals 
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MATHEMATICAL  ABSTRACT 


Several  problems  Involving  estimating  seismic  magnitudes  and 
distance-amplitude  slopes  are  considered  when  there  may  be  severe 
censoring  introduced  by  clipped  or  missed  signals.  It  is  shown  that 
these  classical  situations  as  well  as  the  one  treated  by  Ringdal  (1976) 
are  all  special  cases  of  a  doubly  censored  regression  model.  Maximum 
likelihood  estimators  for  parameters  such  as  station  corrections, 
magnitude,  and  distance-amplitude  corrections  can  be  easily  calculated 
using  a  version  of  the  recently  developed  (Dempster  et  al.,1978)  EM 
algorithm.  Magnitude  estimators  are  compared  for  several  different 
models  and  a  simulation  study  under  conditions  representative  of  papers 
in  the  seismic  literature  indicates  that  biases  of  .2  and  .4  in 
conventional  magnitude  and  slope  (log1Q  amplitude  versus  log10  distance) 
estimators  are  essentially  eliminated  by  the  maximum  likelihood 
procedure . 
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INTRODUCTION 


The  problem  of  estimating  seismic  magnitudes  using  amplitudes  read 
at  a  number  of  recording  stations  is  frequently  complicated  by  the  fact 
that  the  data  may  be  heavily  censored.  This  arises  either  because  of 
clipping  where  all  amplitudes  exceeding  a  given  upper  limit  are  lost,  or 
because  of  a  missed  signal  which  does  not  exceed  a  given  lower  threshold. 
If  one  simply  averages  the  amplitudes  for  those  stations  which  detect, 
without  regard  for  those  stations  not  recording,  serious  biases  may 
result  in  the  estimated  magnitude  levels. 

The  problem  of  estimating  the  mean  and  variance  by  maximum 

likelihood  in  the  presence  of  singly  censored  data  has  been  considered 

in  the  statistical  literature  by  Cohen  (1959),  Swan  (1969),  with  the 

initial  seismic  application  given  by  Ringdal  (1976)  who  exhibited  via 

simulations  the  bias  in  the  conventional  estimator  and  the  reduction  in 

bias  achieved  by  the  maximum  likelihood  estimator  (see  also  von  Seggern 

and  Rivers  (1978)).  Since  the  problem  in  this  form  involves  only  the 

2 

two  parameters  u  and  o  .  simply  scanning  the  likelihood  non-linear 
function  for  a  maximum  in  a  restricted  range  was  used  in  Ringdal  (1976). 

While  the  procedure  described  above  is  acceptable  for  a  limited 
number  of  parameters  lying  in  a  specified  range,  it  cannot  be  used  in 
multiparameter  situations  such  as  those  in  which  station  effects  or 
distance-amplitude  slopes  are  to  be  estimated  simultaneously.  (See,  for 
example,  Richter  (1958),  Nuttli  (1973),  Bollinger  (1973),  Jones  (1977). 
and  Blandford  et  al  (1980)).  Some  early  results  along  this  line  were, 
however,  achieved  by  Christoffersson  et  al.,  (1975).  A  further 
generalization  of  interest  would  be  the  extension  to  doubly  censored 
(clipped  and  missed)  and  completely  missed  observations.  The  large 
number  of  parameters  appearing  in  these  more  detailed  models  make 
classical  non-linear  methods  like  scoring  or  Newton-Raphson  of  little 
use  and  we  employ  a  version  of  the  recently  developed  E-M  ( Expectation- 
Maximizatior;  algorithm  which  was  suggested  by  Dempster  et  al  (1977)  for 
problems  of  this  genre. 

The  general  amplitude-distance  models  that  we  will  use  here  are 
special  cases  of  that  used  by  von  Seggern  (1973)  who  represents  the 
observed  amplitude  at  station  i  for  event  j  as 
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•  •  •  * 


for  i  =  1 


V  j=1 


where  is  the  effect  of  the  ith 


station,  E.  is  the  effect  of  the  ith  event,  and  r.  .  is  the  distance  in 
j  ij 

degrees  from  the  ith  event  to  the  jth  station.  The  parameter  y  gives 
the  rate  of  decay  in  the  amplitude  as  a  function  of  distance  and  is  a 
constant  scale  factor.  It  is  convenient  to  linearize  equation  (1)  by 
taking  logarithms  to  get 


log  =  log  a  +  logS^  +  logEj  “  y  (f^j/10) 


(2) 


which  can  be  converted  to 


'ij 


u  +  ♦ 


e  ♦  6d  +  c .  . 
j  i  j  iJ 


(3) 


where  the  obvious  identifications  are  made  and  e.  .  denote  independent 

1 J  2 

zero-mean  additive  Gaussian  errors  with  common  variance  a  .  The 

"s  „  „  Se 


parameters  are  restricted  to  sum  to  zero,  i.e.  1°  Sj  =  0,  T. 
with  the  magnitude  of  the  jth  event  interpreted  as  the  parameter 


i  -  =  0 


mj  *  v  * 

for  j  =  1,...,n  .  This  is  a  standard  linear  regression  model  of  the 
form  used  in  Appendix  I-A.  The  E-M  algorithm  has  been  used  for  the 
general  linear  regression  model  with  singly  censored  data  in  Aitkin 
(1981)  and  in  Schmee  and  Hahn  (1979);  we  extend  the  results  to  include 
the  possibility  of  double  censoring  and  missing  observations  in  Appendix 
I-A. 

Although  we  will  analyze  several  more  general  versions  of  equation 
(3)  in  the  following  sections,  it  is  interesting  to  note  that  the  E-M 
algorithm  applied  to  the  original  Ringdal  (1976)  model 

ai  1  :  "1  +  E  i  1  (i,) 

for  i  =  1,...,  n  gives  the  (r+1)st  iterate  from  equation  (3)  as 
s 


r  +  1 

U 

1 


r 

w-  1 


(5) 


=  n 


-1  n 


s  r 
T.  w 

i  =  l  il 


where  w 


i  1 


1  if  the  observation  is  present  and 


where 


tlhxL  if  311  -  Cil 

1  r  1 


„r  -< 

il  - 


»(ZJ) 


If  ail  >  C i 2 


»(-z12> 


t  .  . 

U_ 


(6) 


(7) 


The  variance  is  updated  using  equation  ( A 1 6 ) -  Equations  ( 5 )  —  C 7 )  display 
the  updated  estimator  as  an  average  of  observations  when  one  has  them 
and  corrected  values  when  one  only  has  thresholds  (t^  )  and  clipping 
levels  ( ti2)  . 

The  next  two  sections  consider  two  versions  of  equation  (3)  which 
are  important  in  practical  applications.  The  first  is  the  estimation  of 
magnitude  for  a  global  deployment  of  stations,  so  that  B  in  (3)  can  be 
regarded  as  a  known  parameter  (the  B-factor).  The  second  application 
focuses  on  estimating  B  for  a  more  regional  amplitude-distance  study 
(see  Blandford  et  al.  (1981))  without  estimating  the  station  effects. 


Joint  Determination  of  Magnitudes  and  Station  Corrections 

A  global  distribution  of  amplitudes  can  often  be  corrected  for 
distance  using  a  known  value  for  B,  so  that  equation  (3)  might  be 
rewritten  as 


a^j  5  y  ♦  ♦  ®j  ♦  eij»  (8) 

where  the  magnitude  parameter  of  interest  is 

nij  =  y  +  e^  (9) 


1-7 


Table  I  shows  a  typical  set  of  amplitude  for  ng  =  i4  events  which 
theoretically  could  have  been  observed  by  each  of  ng=112  stations.  The 
observations  are  assumed  to  be  in  the  four  categories  described  in  the 
Appendix  with  "<"  denoting  a  category  1  observation  known  only  to  be 
below  the  given  value,  ">"  denoting  a  category  2  observation  known  to 
have  exceeded  the  given  clipping  level  and  "  "  denoting  a  no  show  or 
category  3  observation. 

Again  the  equations  for  the  E-M  algorithm  take  on  a  simple  form  with 


(11) 


estimating  the  station  effect  and  magnitude  value  respectively.  In  this 
case,  the  notation  signifies  that  a  mean  is  to  be  taken  over  the 
appropriate  subscript.  The  values  of  w^  are  determined  by  (A13)  and 
(A8),  where  in  this  case 


(12) 


is  the  estimated  mean  at  the  rth  iterate.  One  can  take  as  initial 
values  the  station  and  event  means. 


Table  II  compares  maximum  likelihood  estimators  computed  under  three 
separate  assumptions,  with  the  original  simple  means  of  the  event 
amplitudes  in  Table  I.  The  means,  shown  in  the  first  column,  are  quite 
high  because  each  missed  signal  puts  a  threshold  value  into  the  average. 
All  distances  0°  <  A <180°  are  used  for  estimation,  and  all  signal  mea¬ 
surements,  even  those  <2  mm,  are  treated  as  signals  and  not  as  noise. 
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TABLE  I 


Observed  P-Wave  Amplitudes  for  Four  Events 


EVENT 


STA 

SHOAL 

PILEDRIVER 

RUBIS 

SAPHIR 

AAE 

<5.44 

<5.72 

5.29 

5.81 

AAM 

5.10 

5.66 

5.29 

_ 

ADE 

<5.99 

<5.96 

<  5.58 

<5.88 

AFI 

<  5.41 

<  5.59 

<5.96 

<  6.14 

AKU 

— 

5.27 

— 

6.08 

ALQ 

— 

5.05 

5.86 

5.97 

ANP 

— 

<6.50 

<7.02 

<7.01 

ANT 

<  4.71 

5.86 

5.44 

5.75 

AQU 

<  4.98 

— 

>5.10 

5.89 

ARE 

5.22 

6.09 

5.45 

5.73 

ATL 

5.12 

5.62 

6.08 

_ 

ATU 

<  6.08 

<6.13 

5.51 

>5.58 

BAG 

<6.14 

<6.03 

<6.11 

<6.10 

BEC 

<  4.93 

<5.10 

— 

_ 

BHP 

<  4.64 

<5.23 

6.01 

6.01 

BKS 

4.99 

5.74 

<  5.79 

<6.09 

BLA 

5.08 

5.50 

— 

5.92 

BOG 

<  4.99 

5.27 

<  5.14 

BOZ 

4.81 

5.36 

5.82 

6.07 

BUL 

— 

5.67 

5.44 

5.58 

CAR 

5.69 

5.85 

5.08 

5.32 

CHG 

_ 

_ 

_ _ 

5.68 

COL 

— 

5.80 

- - - 

4.66 

COP 

<  5.13 

<5.61 

5.44 

<5.56 

COR 

>  5.44 

>5.79 

<6.21 

— - 

CTA 

<  5.94 

<6.11 

5.62 

5.74 

DAL 

— 

— 

5.55 

5.79 

DAV 

— 

<6.99 

— 

<6.47 

DUG 

4.80 

5.34 

CG 

=t 

in 

5.74 

ESK 

— 

5.28 

— 

5.87 

FLO 

4.98 

— 

— 

6.00 

GDH 

<  4.94 

<5.24 

6.01 

6.14 

GEO 

<  5.44 

5.40 

6.08 

5.69 

GIE 

— 

5.57 

— 

GOL 

4.90 

5.56 

5.06 

5.15 

GSC 

5.31 

>4.46 

5.40 

5.70 

GUA 

<  5.56 

— 

<6.39 

<6.39 

HKC 

<  6.42 

<6.60 

<6.18 

_ _  _ 

HLW 

— 

— 

— 

>6.13 

HNR 

<  5.34 

<5.87 

<5.78 

<5.78 

HOW 

— 

<7.17 

— 

_ _ 

1ST 

<5.72 

<5.96 

_ 

_  _  —  _ 

JCT 

— 

5.53 

_ 

, _ —  _ 

JER 

— 

<5.87 

— _ 

6.26 

KEV 

<  5.08 

5.39 

5.68 

6.26 

KIP 

5.47 

6.10 

<6.27 

6.09 

KOD 

<5.63 

1-9 
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TABLE  1  .Continued) 


KON 

4.92 

5.78 

5.38 

5.98 

KTG 

4.88 

5.83 

<4.96 

5.22 

LAH 

<7.41 

<7.25 

— 

<5.64 

LEM 

— 

<5.43 

— 

— 

LON 

4.87 

5.57 

<5.12 

5.63 

LOR 

— 

5.70 

— 

— 

LPA 

<6.56 

<6.09 

<6.27 

LPB 

5 .  24 

— 

5.22 

5.47 

LPS 

- - 

5.41 

5.22 

— 

LUB 

<5.03 

— 

— 

6.45 

MAT 

— 

>5.43 

— 

— 

MAL 

<4.91 

5.35 

6.25 

6.65 

MAN 

<6.20 

<6.72 

<7.98 

<6.88 

MDS 

5.07 

— 

5.38 

5.55 

MNN 

4.42 

— 

5.40 

5.73 

MSH 

— — 

<6.89 

— 

— 

MUN 

<5.61 

<5.90 

— 

5.66 

NAT 

— 

5.84 

— 

— 

NAI 

<5.31 

5.95 

5.43 

>5.86 

NDI 

<5.74 

<5.60 

5.29 

5.64 

NNA 

4.45 

5.64 

5.41 

5-86 

NOR 

— 

5.05 

— 

5.23 

NUR 

<4.83 

5.48 

5.12 

5.42 

OGD 

5.08 

5.  10 

— 

5.66 

OXF 

— 

5.98 

— 

— 

PDA 

<5.86 

<5.38 

— 

<5.93 

PMG 

<5.70 

<5-97 

<5.41 

5.54 

POO 

— 

<5.26 

— 

5-64 

PRE 

_ 

_ 

5.30 

5.43 

PTO 

— 

5.31 

5.32 

5.69 

PEL 

_ — 

5.47 

— 

6.04 

QUE 

<5.16 

<5-88 

4.70 

4.93 

QUI 

<5.71 

<6.00 

<5.93 

<6.23 

RAB 

5.57 

— 

<5.83 

<6.19 

RAR 

— — 

<5.40 

— 

— 

RCD 

5.40 

— 

6.73 

6.45 

RIV 

<6.68 

<6.46 

<5.30 

<5.30 

SBA 

<4.92 

5.27 

5.50 

<5.80 

SDB 

- - 

5.47 

— 

5.78 

SCP 

4.80 

4.97 

5.53 

5.76 

SEO 

4.83 

5.53 

— 

— 

SHA 

<5.15 

5.64 

5.86 

6.16 

SHI 

<5.72 

<5.63 

5.81 

>5.37 

SHK 

— 

5.50 

— 

— 

SHL 

<5.39 

<5.52 

5.33 

5.45 

SJG 

— 

5.86 

— 

5.83 

SNA 

— 

6.06 

— 

— 

SNG 

— 

<5.55 

— 

— 

SOM 

— — 

<6.55 

— 

<6.92 

SPA 

5.11 

5.99 

5.45 

5.82 

STU 

<4.86 

5.47 

5.90 

6.18 

TAB 

— 

<6.49 

— 

— 

TAU 

<5.77 

<6.05 

<5.25 

<5.42 

TOL 

— 

5.55 

— 

— 
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TABLE  II 


Comparison  of  Event  Means  and  Several  Maximum 
Likelihood  Estimators  using  Data  of  TABLE  I 


(1) 

(2) 

SHOAL 

5.24 

4.81 

PILEDRIVER 

5.75 

5.52 

RUBIS 

5.64 

5.46 

SAPHIR 

5.86 

5.75 

a 

(1)  Estimators  computed 
observations . 


(3) 

(4) 

4.84 

4.81 

5.50 

5.51 

5.47 

5.49 

5.74 

5.77 

•  38 

.25 

as  event  means  including  censored 


(2)  Maximum  likelihood  estimators  for  single  event  at  a  time 
using  Ringdal  (1976). 

(3)  Maximum  likelihood  estimators  jointly  for  magnitudes 
excluding  station  corrections. 

(4)  Joint  maximum  likelihood  estimators  for  magnitudes  with 
station  corrections. 
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The  other  columns  are  roughly  comparable  even  though  incorporating 
station  corrections  into  the  computations  reduced  o  (from  .38  to  .25), 
and  the  likelihood  ratio  test  that  the  station  effects  are  0  gave 

x2  =  -2  (-149.82  +  40.17) 

=  219.3 

using  (A17),  where  one  may  compare  with  a  chi-square  random  variable  with 
112  degrees  of  freedom.  Using  the  fact  that  for  large  degrees  of 
freedom  N  the  distribution  is  approximatley  normal  with  mean  N  and  variance  2N, 
we  obtain  2=7  which  is  highly  significant. 

The  iterations  leading  to  the  maximum  likelihood  estimators  in  the 
last  two  columns  take  very  little  time  with  the  100  iterations  required 
in  the  last  case  requiring  less  than  two  minutes  on  the  IBM-360/44. 

Estimation  of  the  Distance-Amplitude  Slope 

For  certain  kinds  of  regional  data  (see  Blandford  et  al.  (1981), 
Nuttli  (1973).  Bollinger  (1973).  Jones  et  al.  (1977)).  The  emphasis 
shifts  to  estimating  a  distance  correction  B  in  (3)  where  one  may  not 
have  enough  data  to  separate  out  a  station  effect.  That  is,  we  consider 
the  model 


aij  =mj  +  Bdij  +  eij  (13) 

2 

for  j=1,...,  ngt  i=1,...t  ni  where  m1 ,  m^ .  mne’  B  and  a  a11  need  to 

be  estimated,  but  we  are  primarily  interested  in  the  effects  of  missed 
and  clipped  signals  on  the  value  of  p.  Again,  the  E-M  algorithm  gives  a 
simple  updating  procedure  of  the  form 


l  (d,,  -  d  .)  (wr  -  wr  1 

ij  ij  • r  1  ij  w-j' 


r+1 


r+1  r 

mj  ■  W.J 


i5  (diJ  " 


”  6r+1d.j 


(14) 
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with  w 


r 


defined  again  by  ( A 1 3 )  and  (A8)  where 


ij 


v 


Vij 


(15) 


Note  that  unless  o  is  fixed,  it  is  updated  using  (A16)  at  the  end  of 
each  iteration.  The  initial  estimator  in  each  case  was  taken  to  be  the 
least  squares  estimated  using  only  the  observed  data  since  this  is  what 
corresponds  to  the  usual  approach  in  the  seismological  literature. 


In  order  to  evaluate  the  possible  improvement  due  to  using  the 
maximum  likelihood  estimators,  consider  the  linear  model  specified  in 
(13)  with  eight  stations  recording  five  events  at  distance  ranges  of 
1 . 3 » 5 .7 , 9. 1 3. 1 7 .  and  21  degrees  respectively.  Suppose  furthermore,  that 
the  upper  and  lower  threshold  are  104  mp  (t2j=4)  and  102  mu  (t1j  =  2) 
respectively,  and  that  amplitude  decays  as  one  over  distance  squared 
(8=  -2).  These  values  are  a  rough  model  of  the  data  seen  in  Nuttli 
(1973)  • 

Table  III  shows  a  collection  of  simulated  amplitudes  for  five  events 
where  the  magnitudes  nij  ,i  =  1 , ..  .,5  were  drawn  randomly  from  a  standard 
magnitude  frequency  law  which  assumes  that  magnitudes  (log1Q)  are 
uniformly  distributed  over  some  arbitrary  interval  (2  to  5  for  this 
example).  (A  constant  of  2  is  inserted  in  equations  14  so  that  both 
magnitude  and  threshold  are  realistic).  The  noise  terms  e^  in  equation 
(3)  were  taken  as  uncorrelated  zero-mean  normal  variables  with  noise 
variance  equal  to  .3  ,  a  typical  value  for  the  standard  deviation  of 
station  magnitudes.  The  theoretical  slope,  as  mentioned  above,  was 
B=-2.  A  further  restriction  was  that  at  least  three  stations  must 
detect  in  order  to  simulate  the  fact  that  an  event  is  not  used  unless 
detected  and  located.  The  effect  of  the  missing  observations  can  be 
noticed  in  Figure  1  as  a  pronounced  tendency  to  flatten  out  (decrease) 
the  negative  slope  of  the  regression  relation  between  amplitude  and 
distance. 


The  estimators  for  the  magnitude  and  slope  parameters  are  shown  in 
Table  IV  where  the  ordinary  least  squares  solution  ignores  the  clipped 
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TABLE  III 


Observed  Amplitudes  at  Eight  Stations  for  Five  Events 
Observed  at  a  Theoretical  Distance-Amplitude  Slope  of  -2.00 


Events 


Station 

1 

2 

3 

lx 

5 

Distance 

1 

5.41 

H 

5.69 

•ft 

5.88 

-1.00 

2 

4.89 

5.04 

4.31 

5.02 

O 

GO 

=T 

-.52 

3 

4.52 

4.79 

4.70 

5.46 

4.65 

-.30 

4 

• 

4.79 

4.63 

5.55 

4.72 

-.15 

5 

» 

• 

4.75 

5.22 

• 

-.05 

6 

• 

4.01 

« 

4.63 

4.03 

.  1 1 

7 

» 

« 

« 

4.24 

• 

.23 

8 

• 

• 

» 

4.50 

4.01 

.32 

• 

Censored 

below,  i.e. 

failed 

to  exceed 

the  lower 

threshold  (4) 

**  Clipped,  i.e.  exceeded  upper  threshold  (6) 
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CUPPING  LEVEL 


lustration  of  effects  of  clipping  and  miss 
servations  to  reduce  estimated  amplitude-d 


and  censored  amplitudes.  The  bias  in  the  estimator  for  the  slope  is 
seen  to  be  substantial;  the  value  is  -1.26  in  this  case  as  compared  to 
the  maximum  likelihood  procedure  which  yields  a  value  of  -2.04. 

TABLE  IV 

Theoretical  and  Estimated  Parameters 


Theoretical 

Least  Squares 

Maximum  Likelihood 

Magnitudes 

m1 

m 

3-54 

4.18 

3.26 

4.20 

4.34 

4.09 

4.16 

4.31 

3.73 

m  3 
m  j. 

4.88 

4.88 

4.84 

m5 

4.13 

4.36 

3-90 

Slope 

6 

-2.00 

-1.26 

-2.04 

In  order 

to  determine 

whether  the  bias 

observed  in  the  above 

particular  case  was  atypical,  an  experiment  consisting  of  drawing  a 
random  set  of  amplitudes  and  performing  seven  iterations  of  the  EM 
algorithm  was  repeated  40  times.  The  bias  terms  in  Table  V  he  low  were 
calculated  by  averaging  the  observed  differences  between  the  true  and 
estimated  values. 

TABLE  V 

Average  Bias  of  Least  Squares  (LSE)  and  Maximum 
Likelihood  (MLE)  Estimators  in  Simulation 

LSE  MLE 

Magnitude  -.20±.02«  .05±.02 

Slope  -.37±.07  -.02 ±.08 

• ±  2  standard  errors 

These  results  indicate  that  the  average  magnitude  and  slope  biases  for 
the  LSE  are  significantly  greater  than  zero  whereas  the  MLE  slope  bias 
is  essentially  zero  and  the  magnitude  bias  is  small.  The  results  for 
the  slope  indicate  that  the  average  of  the  LSE's  should  be  between  -1.70 
and  -1.56  with  95%  confidence  whereas  a  95%  confidence  interval  for  the 


average  slope  produced  by  an  MLE  would  be  -2.08  to  -1.90.  The 
difference  in  slope  is  typical  of  that  seen  as  differences  between  the 
authors  mentioned  in  the  Introduction  for  work  in  the  same  region. 
While  we  do  not  assert  that  any  particular  authors  work  is  in  error  due 
to  the  particular  technique  used,  we  believe  that  use  of  the  MLE 
techniques  will  enable  differences  due  to  geophysics  to  be  more  clearly 
detected . 

We  may  also  compare  the  mean  square  errors  for  the  two  procedures  to 
determine  whether  or  not  the  decrease  in  bias  necessarily  implies  an 
increase  in  uncertainty.  Table  VI  below  gives  the  average  mean  square 
errors  obtained  by  averaging  the  squared  deviations  between  the 
theoretical  and  estimated  values  over  the  40  replications. 


TABLE  VI 

Average  Mean  Square  Error  for  LSE  and  MLE  in  Simulation 


LSE 


MLE 


Magnitude 

Slope 


.087 

.194 


.059 

.081 


Although  it  can  be  seen  that  the  MLE's  have  smaller  values,  the 
reductions  are  almost  entirely  due  to  the  bias  reduction.  For  example, 
since 

p 

Variance  =  Mean  Square  Error  +  (Bias)  , 

we  note  that  the  average  variance  of  the  LSE  of  the  slope  parameter  is 
.055  as  compared  to  the  variance  of  the  MLE  slope  which  is  .081.  Thus, 
it  appears  that  the  variances  of  the  two  estimators  in  this  small  sample 
case  are  roughly  comparable. 

The  small  scale  simulation  given  above  indicates  that  a  rather 
substantial  bias  may  exist  for  conventional  estimators  under  conditions 
which  are  not  unreasonable  in  the  context  of  magnitude  estimation. 
Although  a  broad  range  of  experimental  conditions  have  not  been 
simulated,  it  seems  reasonable  to  infer  that  the  slope  bias  would 
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continue  to  be  present  whenever  there  are  substantial  numbers  of  clipped 
or  missed  signals.  The  simulation  also  indicates  that  the  MLE  tend  to 
reduce  the  bias  in  the  slope  estimator  by  using  the  missing  value 
information  to  steepen  the  negative  slope,  producing  on  the  average,  an 
estimator  with  a  value  closer  to  that  indicated  by  results  from  events 
which  are  recorded  without  excessive  clipping  or  censoring. 


APPENDIX  I-A 

MAXIMUM  LIKELIHOOD  ESTIMATION  FOR  DOUBLY 
CENSORED  REGRESSION  MODELS 


IA-1 


First  of  all,  assume  that  there  are  four  categories  of  observational 
data  which  may  be  available  in  any  given  situation  so  that  the  standard 
regression  model  may  be  written  in  the  form 


i  J 


*ij  1  +  eij '  1  =  °* 1 ,2,3,  j=1 ’ * 


n . 

l 


(AT  ) 


for  the  jth  observation  in  category  i  whereB  is  a  qx 1  regression 
vector,  with  the  qxl  design  vector  _x.j  chosen  to  generate  the 

appropriate  linear  combinations  of  the  parameters  P  ^ ,  P^ . 

The  subscript  i  denotes  the  censoring  category  of  an  observation 
according  to  the  following  class  ications: 


1 . 

i=o  denotes 
available . 

an  observation  where 

Xoj 

and 

yoj 

are 

both 

2. 

i  =  1  denotes 
to  be  below 

that  x  .  is  observed 

the  threshold  t . . . 

1 J 

but 

y1j 

is 

known 

only 

3. 

i=2  denotes  that  x  .  is  observed  but  y 
to  have  exceeded  tn^  clipping  threshold Jt 

1  s 

2  j  * 

known 

only 

4, 

i=3  denotes 

that  x^j  is  observed 

but 

y3  j 

i  s 

not  observed. 

for  notational  ease,  it  will  be  convenient  to  define 


u 


ij 


F 


( A2) 


and  note  that  if  the  errors  e^  are  zero-mean  independent  standard 
normal  variables  with  common  variance  o^,  the  likelihood  functions  can 

be  written,  using  the  fact  the  y . .  is  normally  distributed  with  mean  u 

2  ^  ^  ^  J 

and  variance  a  .  This  means  that  the  likelihood  of  the  incomplete  data 

sampler  yo  =  (  yo  , .  y°  ,no>  .  L,  =  <fcn . -2  = 

(t_v...vt_  can  be  written  as 
Z\  ^ ,nz 


,n  -tl  »-^2  =  _no  2n(27ia2)-- (yoi-UQ.)2 

T~  2o  j=l  (A3) 

nl  n  2 

+  T.  fin  ♦  (Z ,  .  )  +  Y.  In  4>(-Z,,) 
j=l  j=l  J 


1A-2 


7-i1 


(Ay) 


and 

$(u)  =  r  exo{-%  x^}  —  C  A 5 ) 

J  Slv 

—  00 

The  likelihood  function  could  be  treated  in  this  form  by  differentiating 

2 

(A3)  twice  with  respect  to  g  and  a  and  applying  an  iterative 
Newton-Raphson  or  scoring  algorithm  to  estimate  the  parameters.  If 
there  are  a  substantial  number  of  parameters,  this  quickly  becomes 
computationally  intractable,  as  the  corrections  at  each  stage  involve 
inversion  of  the  1/2(q+1)  (q+2)  dimensional  information  matrices. 

A  more  manageable  alternative  procedure  can  be  developed  by 
employing  the  EM  (Expectation-Maximization)  algorithm  as  described  in 
Dempster  et  al .  This  involves  successive  maximizations  of  the 
expectation  of  the  complete  data  likelihood  conditioned  on  the  observed 
data.  For  example,  if  the  complete  data  y^y^.y^.y^,  were  available  in 

the  present  case  with  y.^  =  (yn,yi2 . vi,nl)  for  1  =  the 

complete  data  likelihood  would  be  written  as 

2 

m  l  (y0.y1 .y^.y,  |  °')  = 

2  /  1  2 )  3  ni 

-’-sNln(2no  )(  Jo  IT  Z  (y  .  . 

1  i=o  j=l  1J 

3 

where  u.  .  is  defined  in  equation  (A2)  and  N=  Z  n^  If  the  current 

1 J  2  2  i =  ° 

estimators  for  d  and  o  are  ^  and  o  ,  the  EM  algorithm  defines  ijr+1 

and  .  as  the  values  of  3  anc*  o  maximizing  the  function 

w  r+ 1  - 

Q(B,a2)  =  Er  {Jin  L(  ^  |  g,  02)|  ,  y.,  .  <_t  1  ,j^>  t? }  (A7) 

where  E  denotes  the  expectation  with  respect  to  the  parameter  guesses 
r 

at  the  current  stage. 


I  A- 3 


2 

Now,  by  expanding  (y_  -  u^)  ,  we  note  that  all  the  terms  in  the 
conditional  expectation  given  in  (A6)  can  be  determined  when  we  know 


Er  O'ljK,  -  ■  hj 


♦<2lt> 

‘<2Ij> 


(A8) 


and 


Er  'hjlhj  >  hj>  ■  hi  +  ”, 


*(Z21> 

t(Zl i> 


Er  'hlihi  <  hi>  -  hj2  + 


♦<ZU> 


(A9) 


Er  ‘hi'lhj  >  hj>  -  “2i  «r+h  (<=2j  +  l”Jj> 


where  u r  and  Zr  are  equations  (A2)  and  (A4)  with  g  and  0  2 

ij  2  -r  r 

replacing  p  and  a  .  The  function  $  (•)  is  defined  as  the  standard 

normal  density 


#  (u)  =  (?*)  i  e*P  i-  h  u?} 


(  A 1  0 ) 


Then,  returning  to  (A6)  and  (A7) t  we  note  that  the  succeeding 
estimators  can  be  determined  by  maximizing 


■i,2  J-l  Erf<hj  -  hi’2!2!,  -  hi” 


Ja 


"lo2  J=1  Er{(y2j-U2j)2|y2j>t2j)} 


'2a2  J=1  Er{(y3j~P3j)2} 


(All) 


IA-4 


2 

with  respect  to  6  and  a  . 


Now  in  order  to  maximize  the  above  expression,  note  that 
written  in  the  form 


Q(.5,o2)  =  \  N  ln  (2na  ?) 


—  ,,[£  vr.  -  2  6'  E  x.  w^.+fe'  7. 


-  l  —  V  .  . 

2a2  ij  U 


ii-iJ  ~ 


i  J 


x  x  6  ] 


where 


r 

w  .  . 
ij 


and 


i  =  0 


Er(^ij  y i j—1 ij  3 ’ 


ij  - 


yij ’  1-° 

Er<yijly£ji‘ij>>  ‘-1 
Er<yijlyij>tij>'  *-2 

i^V^v  1'3 


with  the  conditional  expectations  given  in  (A8)  and  (A9).  Some 
2 

shows  that  Q(§,o  )  is  maximized  for 


and 


— r+! 


=  2i4  *  ^ )  _1  E  x  .  wr  . 

ij  ^  ij_iJ  ^ 


r  +  1 


n-1  rno 

N  2  (y 

j  =  l 


Oj 


r  \  2  ,  2 

‘l,0j)  +° 


r  [N"V^1!lf2i)Zr 


'♦(Z2J) 


2j 


E1  ^(Zr. 

j  =  1?Tz^ 


1 


1A-5 


it  may  be 


( A12) 


(  A  1  3  ) 


(  A  1 4  ) 


algebra 


( A15) 


(  A 1 6 ) 


It  is  convenient  to  notice  that  ( A 1 5 )  is  just  the  ordinary  least  squares 

or  maximum  likelihood  solution  evaluated  at  the  observations  wr .  .  as 

ij 

defined  in  ( A 1 3 ) -  This  facilitates  the  application  of  the  technique  to 
the  special  cases  described  in  the  text. 

The  basic  iterative  procedure  can  be  started  with  initial  estimators 
2 

for  8  and  o  derived  from  maximum  likelihood  estimators  derived  by 
fixing  the  censored  values  at  the  upper  or  lower  thresholds.  Then, 
equations  (A15)  and  ( A 1 6 )  can  be  applied  to  generate  successive  iterates 
for  8^  and  o  until  the  incomplete  data  log  likelihood  given  by  (A3)  is 
maximized . 

It  should  be  noted  that,  for  reasonably  large  samples  certain 
hypotheses  can  be  tested  by  comparing  the  log  likelihoods  computed  under 
the  various  models.  For  example,  if  estimates  n ^  parameters  and 
is  a  restriction  of  Hq  which  estimates  n^  parameters  then 

v2  =  -2 (In  L  -  In  L  )  (A-17) 

*  u  n 

is  distributed  as  chi-square  random  variable  with  (n^-n^)  degrees  of 
freedom  if  H  is  true. 


APPENDIX  II 

CRUSTAL  MODELS 

« 


II-l 


Station 


AAE 


2.7  Fair 


3.6  2.7  i  Poor 


Good 


3.3 


3  6 


B.11  3J) 


2.0 


1.0  2.0 


0.8  3.0 


Crustal  Models  _ 


Geology 


1000  m.  Tertiary  basalts 
on  1200  m.  Mesozoic  sandstone 

and.  U®  £  si  Qri.e_Qn_b_a£  e  rfi  e  ijt _ 

glacial  drift 
shale 

limestone,  dolomite 


PreCambrian  deposits 


basalt 


PreCambrian  granite 


Pleistocene  andesite 

.  laya.fi  _ 

limestone 


limestone  over 


Reference 


U.  of  Mich. (1964) 


U.  of  Mich.(1964) 


U.  of  Mich. (1964) 


U.  of  Mich(l964) 


Noponen  and  Cass  (1930) 


U.  of  Mich. (1964) 


U.  of  Mich.  (1964) 


U.  of  Mich. (1964) 


U.  of  Mich.(1964) 


80  m.  unweathered  volcanics 
fivey  b55.exa.ei1l _ 

U.  of  Mich.  (1964) 

PreCambrian  migmatites 

King  and  Beikman(1974) 

1.0  6.1 


Mesozoic  rocks? 


limestone 


aeolian  limestone 


(no  info.,  used  nearby 
model) 


velocity  model 


dolomite 

shale 

dolomite 

.  PK-Csmbosa  s^diiBSois _ 

sandstone 


velocity  model 


lava,  used  nearby 

.msdsi - 

Cretaceous  schists 

.fitflWiS - 

granite 


PreCambrian  sediments 

.  SUd-XfilSAIUPS _ 

PreCambrian  metamorphics 


Choubert  and  Faure-Muret(1980)  : 


U.  of  Mich. (1964) 


U.  of  Mich.(1964) 


Steinhart  and  Meyer(1961) 


Steinhart  and  Meyer(l961) 


Lowry(1971) 


U.  of  Mich.(1964)  j 


Steinhart  and  Meyer(1961) 


Steinhart  and  Meverfl961), 

.y^9fJ4i£h.(Ji64jl__ _ 

U.  of  Mich  (1964) 


King  (1969) 


King(1969) 


U.  of  Mich.  (1964) 


IP 


Poor 


Tertiary  (limestone?) 


Station 


COR 


HN-ME 


HNR 


6_1  i  3_6 


6_1_  ■  3_6_  !  2  7 


1.0 


1.0  4.0 


1.0 

«  Le,Ll-3J_L2.7 


JLL13.fi 


_  Crustal  Models 


Quality  Geology 


Good  weathered  basalt, 
velocity  model 


granodiorite 


Reference 


Berg  et  el  (1966), 
U  of  Mich. (1964) 


U.  of  Michigan(1964) 


Good  4150' Cretaceous  sediments  U  .  of  Mich.  (1964) 

- 2iL  Sale  _ 

Good  Paleozoic  basalt  Choubert  and  Faure-Muret(l9B0) 

Fair  granite  *  U  of  Mich. (1964) 


Choubert  and  raure-Muret(1990) 


Permian  continental 


clay  on  Mississippian 
limestone,  velocity 


PreCambrian  diorite 


1.9  j  Poor  |  volcanic  tuff 


granite  and  diorite 


|  Cretacous  limestone 


2.7  !  Good  dolomite  and  slate 

_ .granite _ 

Poor  sandstone  on 

_ Mipc  en_e_limytqne  _ 

Poor  alluvium 


Devonian  graywacke 


2.6  |  Poor  lower  Cretaceous 


Poor  Upper  Cretaceous 


granulite 


Howe  and  Koenig(  1 961 ), 

Stauder  et  al(l901), 

_  U  _of  _Mi  ch.  (J  964J _ 

U.  of  Mich. (1964) 


U.  of  Mich. (1964) 


Choubert  and  Faure-Muret(lS30) 


Noponen  and  Cass(1980) 


ECAFE(1971) 


U.  of  Mich.(1964) 


Teledyne  Geotech(1966) 
U.  of  Mich. (1964) 
j  Noponen  and  Cass(190O) 


U.  of  Mich.(1964) 


King  and  Beikman(l974) 


Choubert  and  Faure-Muret(  1980) 


U.  of  Mich.(1964) 


Poor  alluvium 


granite 


gneiss 


Noponen  and  Cass(l280) 


Holtedahl  and  Dons(1960) 


Noponen  and  Cass(1930) 


33. 

.32.. 

, _  - _  ! 

3.6 

2.7 

Fair 

Tertiary  basalt 

U.  of  Mich. (1964) 

I - T 

IStation  p 


j _ JtasJsL. 

it_J.jl_L.0_ 


|LAH  [  0.3  Ti.oT l.o"~2.0  "  Fair" 
_ _ 


Quality 


I _ _  i.  2 ■+ _ 

LEM  1.0  [ 6.4  3.2  2.8  Poor 

-  6.4  3.2  2.8 


6  1 

3.6~ 

6.1 

3.6 

1.0  !  2.0 


1.0  |  3.6  2.0  2.5  Poor 


1.0  6.1 


|  1.0 

3.6 

-S-iJ 

2.0 
_3.6  J 

2.4 
_2_7_  J 

Poor 

!  i.o 

5.5  ! 

3.2  1 

2.6  ! 

Poor 

10  6.1 


3.6  2.7 


2.4  Poor 


6J_L3.6 


.  Crustal  _Models_ 
Geology 


1200'  alluvium,  bedrock  is 

probably  schist _ _ 

volcanies? 


Middle  Jurassic’’ 


600'  Quaternary  loess 


clay  and  sand 


1 

3.6 

2.7  | 

Fair 

1 

granite 

Pliocene  continental 

deposits _ 

Tertiary  limestone 


Tertiary  tuff  over 
JB&tJcX9i£dQi£  S.  gn.d_.bas  erjj  gnt_ 
Quaternary  extrusives? 


quartzite 


lavas  and  pyroclastics 

on  basement  gneiss _ 

Paleocene  and  Neocene 


quartzite  and  sandstone 


Pennsylvanian  to 

Cretaceous _ 

PreCambrian  gneiss 


ash  and  scoria 


Cretaceous 


Reference 
U.  of  Mich.(1964) 

ECAPE0971) 


U.  of  Mich.(1964)  ! 

i 


Choubert  and  Faure-Muret(1930) 


U.  of  Mich.(1964) 


of  Mich.(1964) 


U.  of  Mich.(1964) 


King  and  Beikman(1974) 
U.  of  Mich.(l964) 


U.  of  Mich. (1964) 

Choubert  and  Faure-Muret(1960) 


U.  of  Mich.  (1964) 


of  Mich. (1964) 


ECAFE(1971)  | 


Noponen  and  Cass(1980) 
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